Universality in the pair contact process with diffusion 
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The pair contact process with diffusion is studied by means of multispin Monte Carlo simulations 
and density matrix renormalization group calculations. Effective critical exponents are found to 
behave nonmonotonically as functions of time or of system length and extrapolate asymptotically 
towards values consistent with the directed percolation universality class. We argue that an inter- 
mediate regime exists where the effective critical dynamics resembles that of a parity conserving 
process. 
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I. INTRODUCTION 

Out of equilibrium systems may display phase tran- 
sitions analogous to those found in their equilibrium 
counterparts. These transitions are classified into dis- 
tinct universality classes, characterized by a set of criti- 
cal exponents [ll Q • Particular attention has been paid 
to one-dimensional systems with transitions from an ac- 
tive state into absorbing states, i.e., frozen configurations 
from which the system cannot escape. For these systems, 
so far, only two distinct universality classes have been 
firmly established: the directed percolation (DP) and the 
parity conserving (PC) universality class. While the for- 
mer is ubiquitous and found in several systems with very 
different dynamical rules, the latter is only known to oc- 
cur when extra symmetries are present . 

A model which has attracted quite some interest re- 
cently, because it may indeed belong to a novel univer- 
sality class is the so-called pair contact process with 
diffusion (PCPD). Despite a rather intense activity in the 
past couple of years HESSIllillSlillll, the under- 
standing of the PCPD and its relation with other known 
models is still unsatisfactory. In the fermionic version of 
the model — the one studied here — each lattice site is ei- 
ther occupied by a single particle {A) or empty (0). The 
reactions are 



AAO AAA 
OAA AAA 
AA 00 
AO ^ OA 



with rate 



(i-p)(i- 



with rate p 
with rate 



(l-d), 



(1) 



with 0<d<l,0<p<l. The analysis of the critical 
properties of the PCPD has shown to be much more dif- 
ficult than all similar models analyzed so far and several 
scenarios have been proposed. First, a similarity of the 
exponents P/h'± and z — v\\/v±^ (where and are the 
correlation length exponents along the time and space di- 
rections and (3 the order parameter exponent) with those 
of the PC class 4j was reported, although in the PCPD 
there is no conservation of parity. It was later suggested 
that the PCPD could belong to a new universality class 
with exponents close to, but different from the PC val- 
ues 0. It has also been argued that there could be two 



universality classes 1^ at small and large diffusion rates 
(d), or continuously varying exponents J^, or that scaling 
may even be violated [g. More recently the option of a 
slow crossover to DP was also discussed although the 
general belief is that the PCPD belongs to a novel univer- 
sality class [1 H ES m 113 • 

Field theoretical methods, 
which have been successfully applied to other reaction- 
diffusion models failed so far to clarify the critical 
properties of the PCPD Hill- 

In this paper we present some insights into the PCPD. 
We show that accurate numerical results from Monte 
Carlo (MC) simulations and density matrix renormal- 
ization group (DMRG) calculations convincingly demon- 
strate that for sufficiently long times and system lengths 
the exponents crossover towards the DP values. Correc- 
tions to scaling are, however, rather strong and only an 
accurate extrapolation of the effective critical exponents 
allows to identify the final asymptotic critical behavior. 
We give evidence that in the near-asymptotic region the 
model shows effective exponents close to the PC class 
values, which suggests that the critical behavior of the 
system is described by two competing fixed points. 



II. RESULTS ON BULK PARTICLES AND PAIR 
DENSITIES 

Our MC simulations exploit a technique known as mul- 
tispin coding |15| |. The basic idea is that in a simulation 
of 64 systems, each with L sites, the occupation of site 
i in the fcth simulation is stored in the fcth bit of 64-bit 
word A[i\. To perform the reaction AAO AAA in all 
64 systems at a randomly chosen site i and its neighbors, 
one logical operation 1] = 1] V {A[i] ^A[^- 1]) 
suffices, where V and A are the logical operations OR and 
AND, respectively. Other reactions require only slightly 
more elaborate logical operations. A direct implementa- 
tion along these lines might result in 64 simulations, each 
of which statistically correct, but strongly correlated to 
each other because the site selection is shared. To al- 
leviate this correlation without sacrificing efficiency, we 
employ random bit patterns that decide which reaction 
will be attempted in which system. The strong point of 
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FIG. 1: Ratio of the particle density p and pair density p* , as 
a function of t-^' , with 7' = 0.33, for the PCPD with d = 0.5 
and at the estimated critical point p = 0.152 45. This ratio 
approaches a constant. 



multispin coding is its efRciency. As illustrated above, 
only a few logical operations (each usually carried out in 
a single clock cycle without delay) suffice to update a site 
in 64 systems simultaneously. For each combination of d 
and p we simulated 64 systems with i = 100 000 sites 
over 3 x 10^ Monte Carlo time units, in about 15 hours 
on a single-processor workstation. We also simulated 
16 X 64 systems with L = 100 000 sites over 10'' Monte 
Carlo time units on a parallel computer, for d — 0.5 and 
p = 0.152 45, our estimate for the critical point. 

A standard procedure to obtain critical exponents by 
means of Monte Carlo simulations is to study the de- 
cay of the particle density p{t) as a function of time 
starting from a random configuration of particles. At 
the critical point one has p{t) ^ t^P/'^w . To monitor 
the decay it is convenient to define the effective expo- 
nent (5eff = —d\np{t)/d\\it. Typically JeS is plotted as a 
function of 1/t. At the critical point, in the limit i — > 00, 
it approaches a finite value (Jeff ^ while it de- 

viates upwards or downwards with respect to this value 
in the inactive and active phases, respectively. This cri- 
terion allows to estimate the critical point location and 
the ratio P/v\\ [11 • Some care, however, has to be taken 
when corrections to scaling are particularly strong, for 
example when p{t) ~ (1 + ct^f) with 7 < 1 and 

c a constant. In this case 5eS plotted as a function of 
1/t approaches P/fw with an infinite slope. The ideal 
situation would be to plot S^s as a function of 1/V as 
in that case the approach to the asymptotic value would 
be linear. Further on in this paper we will give numeri- 
cal evidence that the correction-to-scaling exponent 7 is 
close to P/vw (which is much smaller than 1). A natural 
choice is therefore to plot 5cS as a function of the particle 
density p instead of 1/t, to avoid infinite slopes. 

Besides the decay of the particle density p, we also 



FIG. 2: Plot of (5cff at d = 0.5 and p = 0.1524, p = 0.152 45, 
and p = 0.1525 (from bottom to top) obtained from the de- 
cay of the particles and pairs densities. The dashed lines are 
parabolic fits to the data. The two horizontal lines show the 
ratio P/v\\ for the DP and PC class. Inset: d = 0.1 and 
p = 0.1110, p = 0.111 05, p = 0.1111. 



consider the decay of the pair density p* = {AA) . Before 
presenting the results for the critical exponents we ana- 
lyze the behavior of the ratio p/ p* at the critical point. 
Such a quantity, which is shown in Fig. ^ for d = 0.5, 
approaches a constant value (sa 2.45) asymptotically for 
long times. The most important consequence of this fact 
is that in PCPD one can extract the critical exponent 
P/i'W both from the decay of the particle and pair densi- 
ties Numerically the ratio is much better behaved 
than the individual densities p and p*, as fluctuations in 
the individual densities are highly correlated and largely 
cancel each other in the ratio. 

The data of Fig. also provide an estimate of the 
leading correction term in the asymptotic limit t — ^ 00, 
which appears to be of the type p/ p* ~ C{l+Dt~'^ ) with 
7' ~ 0.33, as shown by the linear approach to the asymp- 
totic behavior of the data when plotted versus t~'^-^'^. 
This exponent is much larger than its equivalent 7 in the 
particle and pair densities, which implies that the lead- 
ing corrections for p and p* cancel in the ratio p/ p* . We 
tested that a similar cancellation also occurs in the pro- 
cess A — > 2>A, 2A 0, which belongs to the PC univer- 
sality class IJ: the leading correction in p and p* scales 
as t~"'^, whereas the leading correction in the ratio p/p* 
scales as 1/t. 

Figure 121 shows a plot of Scff at d — 0.5 as a function 
of p and p*, calculated in the PCPD from the decay of 
particles and pairs, respectively. Three different values 
for p have been plotted around the critical point which 
we estimate as Pc = 0.152 45(5). For p > pc (inactive 
phase) and p < pc (active phase) Scs rapidly veers up and 
down, as expected. At the critical point. Sag approaches 
the y axis with a finite slope, indicating that the leading 
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FIG. 3: Log-log plots of p (particle density) vs t. (a) d — 
0.5, p = 0.1524, p = 0.152 45 and p = 0.1525 (from top 
to bottom), (b) d = 0.9 and p = 0.2330, p = 0.2335 and 
p = 0.234 (from top to bottom). The critical point densities 
are plotted as thick solid lines. The dashed lines are linear 
fits to the data, shifted for clarity. 



correction to scaling is most likely described by an ex- 
ponent roughly equal to S itself. A parabolic fit through 
the data yields as a common estimate for the critical ex- 
ponent = 0.17. Similar calculations were repeated 
for other values of the diffusion coefhcient d = 0.05, 0.1, 
0.2, and 0.9 (the inset of Fig. |2]shows the case d = 0.1), 
with the same results which we can summarize with the 
estimate P/i^w — 0.17(1). This value is consistent with 
the DP class exponent f3/iy\\ — 0.159. 

An alternative way to calculate critical exponents 
would be to use a linear fit of the densities versus time on 
a double-logarithmic scale. However, such fits are haz- 
ardous when dealing with very slow convergence, as is the 
case here, and may lead to wrong estimates for the crit- 
ical exponents. We illustrate this in Fig. O which shows 
the plot of the average particle density as a function of 
the time in a double-logarithmic scale for d — 0.5 and 
d — 0.9. In the former case a straight line (dashed) fits 
extremely well the critical density decay leading to the 
estimate i5cff — 0.219. The analysis of the effective expo- 
nent (see Fig. |2Jl, however, provides a closer inspection 
of the local slopes of the double-logarithmic data. This 
analysis reveals some remaining curvature, and the final 
estimate of the exponent is significantly lower, compared 
to that obtained from the fit in the double-logarithmic 
scale. In the case of higher diffusion [see Fig. |3Ib)] 
the curvature is more pronounced and clearly visible also 
in the double-logarithmic plot, which can be fitted by 
two straight lines with slopes Scs ~ 0.277 in the range 
4 < Ini < 8, and with ScS « 0.212 for 10 < Int < 15. 
Notice that the former exponent is consistent with that 
expected for the PC class {6pc = 0.286 Q). 

Again an extrapolation of the effective exponent (as 
done in Fig. [21 shows convergence to a value consis- 
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FIG. 4: Plots of 7eff for d — 0.5 at the system critical point 
{p = 0.152 45) for lattices up to L = 60 calculated from 
the decay of particle (circles) and pair (squares) densities as 
function of the lattice lengths. The dotted and dashed lines 
are fits in powers of the densities. Inset: 7eff for d — 0.2. 

tent with DP. Note that the value /3/i^|[ ~ 0.21 is con- 
sistent with the most recent Monte Carlo estimates for 
the PCPD d m 113. In particular, Kockelkoren and 
Chate, ^3 performed a series of Monte Carlo simula- 
tions for a bosonic version of the PCPD where the con- 
straint of one particle per site is released. Their estima- 
tion of critical exponents is based on a straight-line fit 
to a double-logarithmic plot of p versus t, from which 
they find = 0.200(5). This bosonic version of the 

model is claimed to suffer less from corrections to scaling 
than the fcrmionic case. Notice, however, that also in 
the fermionic PCPD studied here the density decay at 
d = 0.5 (see Fig. |2Ia)) is rather straight in a double- 
logarithmic plot for simulation times similar to those in 
Ref. JJJ . The advantage of the effective exponent analy- 
sis performed here is that it allows to extrapolate the nu- 
merical results to time scales beyond those actually sim- 
ulated. 

Next we present some DMRG results. DMRG [l8| al- 
lows to calculate accurate stationary state probabilities 
for chains of moderate lengths As usual in DMRG, 
we used open boundary conditions. In the PCPD on 
a lattice of finite length there are only two stationary 
states: a state with no particles and a state occupied by 
a single diffusing particle. To induce a finite density of 
particles we added a reaction — s- A at the two boundary 
sites. The particle density decays from the two bound- 
aries and forms a U-shaped profile. For chains of various 
lengths we calculated the density of particles p{L) and of 
pairs p*{L) at the central site of a system of length L. 
At the critical point these quantities decay in the limit 
L ^ 00 as p{L) - p*{L) - L-f^/"^. 

Figure ^ shows the effective exponent jcS = 
—d\np{L)/dlnL versus p for d = 0.5 at the critical 
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FIG. 5: Plot of the effective exponent Zeff as a function of 
1/L for d = 0.2, 0.5, 0.8, and 0.9. Inset: blowup of the data 
for d = 0.5. 

point. As in Fig. |21 we include also the data for the 
pairs. Dotted and dashed lines are fits with polynomi- 
als in the densities. Again, a test of good convergence is 
that both exponents extrapolate to the same asymptotic 
value. This requirement seems indeed to be fulfilled and 
we find as extrapolation [3/v± = 0.27(4). This exponent 
is again consistent with the DP value = 0.252 0. 

Similar results have also been found for other values of 
the diffusion coefficient d. Extrapolations for d = 0.2 are 
shown in the inset of Fig. 21 At small d the maxima in 
7efi shift to longer L, thus extrapolations are somewhat 
less stable. In this case we take the estimate obtained 
from the particles f3/i^± = 0.28(5). 

Previous DMRG results j| were restricted to the den- 
sity of particles and to smaller systems than that studied 
here. The data for the effective exponent showed a mono- 
tonic behavior (except at very strong diffusion) and were 
analyzed using an extrapolation with polynomes in 1/L. 
These extrapolations lead to a value consistent with the 
PC class exponent = 0.50 Q. The present cal- 

culation, extended to the density of pairs and to longer 
systems, reveals that nonmonotonicity in the effective ex- 
ponent is a common feature at all d. This nonmonotonic- 
ity leads to a rather strong decrease of the extrapolated 
exponent compared to the estimates of Ref. 0| . 



III. RESULTS ON THE DYNAMICAL 
EXPONENT 

From the ratio of the exponents /3/i'± and P/i^w one 
can estimate the dynamical exponent z = v^i/u^. Since 
both (3/i'± and P/i'w are consistent with DP, also the dy- 
namical exponent z agrees with the DP value zdp = 1.58. 
It is however instructive to show the results of an inde- 
pendent calculation of z. This quantity can be obtained 
from a finite size scaling analysis of A, the gap of the 



FIG. 6: Scaled particle density for d — 0.9 at the estimated 
critical point p = 0.2335 for L = 100, 200, 500, 1000, 2000, 
and 5000. Inset: plot of Inr vs InL for two values of the 
diffusion constant. We estimate z = 1.61(3) for d — 0.9 and 
2 = 1.70(3) for d = 0.2. 



Master operator, which is the inverse of the relaxation 
time of the system (see Ref. 19] for details). As a func- 
tion of the system length L the gap decays as A L~^. 

Figure \S[ shows a plot of the effective exponent z^s = 
— 9 In A/9 In i versus l/L. The calculations are similar 
to those reported in Ref. 0|, but now for longer sys- 
tems (up to L = 46 compared to L = 30 of Ref. jj]). 
The critical point locations were obtained from Monte 
Carlo simulations, which for this purpose are faster and 
more efHcient than DMRG. Therefore we concentrated 
our computational efforts on a single value of p — Pc and 
could obtain results for longer systems. As is clear from 
Fig. |5l the exponent Zoff is rather sensitive to the value of 
the diffusion rate d. As the estimates of /3/j^± and /3/i^\\ 
are instead rather stable as a function of d we contribute 
this sensitivity to rather strong finite-size effects. Notice 
that the finite L corrections change sign from the weak to 
the strong diffusion regime. The border value is around 
d = 0.5 where Zoff has a very weak dependence on L. 
The data (see inset) run extremely close to the PC value 

zpc = 1.75. At higher diffusion rates d « 0.8 0.9 

the effective exponent Zcff for the range of sizes investi- 
gated is much lower than zpc- At the strongest diffusion 
investigated, extrapolations with different forms for the 
correction to scaling terms as 1/L or 1/ VL yield values 
in the range 1.5 ^ z ^ 1.65, which should be compared 
with the DP value zdp = 1.58. Current Monte Carlo 
estimates from various authors 

HIE EH place the expo- 
nent z in the range 1.7 — 1.8 and the calculations were 
mostly performed in the weak diffusion regime. 

We also performed a series of Monte Carlo simulations 
to calculate the exponent z using finite-size scaling anal- 
ysis. At the critical point and on a finite system the 
particle density decays as p = /(tL~^), with / a 

scaling function. For finite L, p follows a power-law de- 
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FIG. 7: Plots of 7^ff for three different values of the diffusion 
constant. Iforizontal lines are the reference exponents for the 
DP and PC universality classes. 



FIG. 8: Plots of vs p" for different values of the diffusion 
constant and absorbing boundary conditions. Inset: 7|f{ vs 
p*^ calculated from the surface pair density. 



cay up to a characteristic time r after which it drops 
exponentially. One expects that r scales as r ~ L^. To 
estimate r we calculated p for lattices up to L = 5000 
and t = 10^ Monte Carlo time units. Figure shows a 
plot of \n.{t^/'^\\ p) versus Ini for various L and d — 0.9. 
The intersection of the data with a horizontal line at 
tl^/'^w p = k (with k a constant) provides an estimate of 
T. As we work in a region where the particle density 
is rather low and fluctuations are large, and as the cal- 
culation of z requires very smooth data, we performed 
averages over a large number of samples (> 10"^). For 
the calculation we used P/v\\ = 0.17, which is the value 
determined above, and k = —2 (see Fig. E)). The inset 
shows a double- logarithmic plot of r versus L for d = 0.2 
and d = 0.9 at their critical points. In the former case 
we restricted ourselves to L = 2000 as the relevant times 
are typically longer at weak than at stronger diffusion, 
as expected. Notice that in both cases the data are well 
fitted by straight lines yielding the estimates z = 1.70(3) 
for d = 0.2 and z — 1.61(3) for d = 0.9, where the lat- 
ter value is consistent with the dynamical exponent of 
directed percolation zup = 1-58. The results generally 
confirm the DMRG findings according to which the dy- 
namical exponent is generically smaller, for finite L, at 
higher diffusivity. We also notice that by varying the 
value of entering in the y axis of Fig. |H|one changes 
the estimate for z. For instance, if we take f3/iy\\ = 0.20, 
as calculated in Ref. this leads to an increase of 0.03 
in the estimated value for z. The estimate of z is rather 
stable for changes in the constant k. 



IV. RESULTS ON SURFACE DENSITIES 

Boundary quantities are easily accessible in DMRG 
techniques |22[, as one is basically forced to work with 



open boundary conditions. Surface criticality in absorb- 
ing phase transitions has been the subject of several stud- 
ies in the past years both for models in the DP j23| and 
in the PC 22] universality classes. In the latter case, it is 
known that there are two distinct surface exponents de- 
pending on the type of boundary conditions applied . 
The results of a DMRG calculation of the surface crit- 
ical exponents for a reaction-diffusion model in the PC 
class are presented in the Appendix. Here, we report on 
the surface critical exponent calculations for the PCPD, 
using the same types of boundary conditions as in the 
Appendix. 

As in the calculation of the bulk particle density of the 
preceding section we inject particles through the reac- 
tion — + A at the boundary site labeled by the position 
i = 1 in order to induce a finite density of particles in 
the system, and we measure the particle density p^{L) 
at the opposite boundary site i = L. Asymptotically 
for L — > oo, we expect p^{L) ~ L~'^°/'^^, where (3^^ is 
the order parameter surface exponent. The two different 
boundary conditions (BCs) applied at the site i = L are: 
(a) No particles are allowed to leave the system from the 
boundary site and (b) particles may diffuse out of the 
system, i.e., the reaction ^ ^ (with rate d) is added 
at that site. We refer to these as reflecting and absorbing 
boundary conditions, respectively. 

In Fig. [71 we plot the effective exponent 7^^^ = 
—d\np^{L)/d\nL versus p" in the case of reflecting BCs. 
Horizontal lines show the ratio (3^ /vi_ for DP (= 0.667 
[Hill) and PC (= 0.72 ^). In the DP case the dif- 
ferent BCs produce the same critical exponent. Effective 
exponents in this case grow monotonically, contrary to 
what is found for bulk exponents. Notice that a cubic 
flt yields a quite stable estimate [3'^ /v±_ = 0.72(1) in the 
range d ^ 0.5, a value actually consistent with the sur- 
face exponent for the PC class (see Appendix). Only 
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at higher d we observe some deviation from PC. The 
fact that the extrapolated surface exponents vary with 
d, while our current estimates for the bulk exponents are 
independent on d, is an indication that the former are 
not yet the correct asymptotic ones. 

Figure |H1 shows 7^^^ versus p'* for the case of absorbing 
boundary conditions. Again for weak diffusion the expo- 
nent seems to extrapolate rather convincingly to values 
close to the PC class {(i'^ /v±_ ~ 1.11, see Appendix) , while 
for strong diffusion it increases to much larger values. 
Also in this case there is no clear signature of nonmono- 
tonic behavior, except for the case d — 0.5 where the 
data for the largest systems pass through a maximum. 

We also analyzed the effective exponent data from the 
pair density p*" which are shown in the inset of Fig. |S1 
in the case of absorbing BCs, and plotted as functions of 
p*^ . In the range d ^ 0.5 the data extrapolate close to 
the PC surface exponent (3'^ /v±^ k, 1.11 as for the parti- 
cle density. At very strong diffusion {d = 0.9) the sur- 
face effective exponent shows a nonmonotonic behavior 
with a maximum around ~ 2.2. Notice that parti- 
cle and pair exponents in this case are rather far apart 
from each other and it is quite hard to find a common 
extrapolation value. We would expect for a similar 
behavior as for the bulk exponents, i.e., an increase fol- 
lowed by a decrease towards the asymptotic value. We 
suspect that in the present surface exponent calculation 
the decreasing side has barely been reached. So we tend 
to distrust the extrapolation as estimates of the genuine 
asymptotic behavior. They rather provide some insight 
on the preasymptotic region and actually point to a sim- 
ilarity with PC surface exponents at weak diffusion. 



V. DISCUSSION 

To conclude, by combining Monte Carlo and DMRG 
calculations we analyzed the critical properties of the pair 
contact process with diffusion. This model has been the 
subject of increasing attention in recent years. Although 
the debate around it has not yet been settled, the main 
belief is that the PCPD belongs to a novel universality 
class which differs from the known DP and PC classes. 

In our opinion, however, the most plausible scenario 
for the PCPD is that it ultimately falls into the DP 
universality class. The asymptotic behavior is, however, 
masked by rather strong finite size and time effects, char- 
acterized by small correction-to-scaling exponents, as our 
Monte Carlo simulations for the decays of the particle 
density p, the pair density p* , and the ratio p/p* have 
demonstrated. 

The exponents and f3/i'± extrapolated both from 
p and p* appear to be stable as functions of the diffusion 
constant d and actually consistent with the DP class val- 
ues. The data show a nonmonotonic behavior both in 
time and system size, which in our opinion points to a 
crossover phenomenon between two competing types of 
critical behavior. The surface exponents, which we also 



investigated, turned out to be instead rather sensitive to 
the value of d; a sign, in our opinion, that the extrapo- 
lated values are probably not the true asymptotic ones. 
Interestingly enough, particularly at weak diffusion, the 
extrapolated values are rather stable and consistent with 
those expected for the PC class. 

In early numerical studies of the PCPD IM IS > re- 
stricted to shorter simulation times and system lengths 
compared to those considered here, several quantities as 
and h'\\/i'± were found to be quite consis- 
tent with the PC class values. It is now generally agreed 
that the PCPD does not belong to the PC universality 
class, as more extensive simulations performed by several 
groups have shown convincing 

ly Hiiliaia stm 

one would like to understand if the observed similarity 
with the PC exponent is purely fortuitous or if there is 
some deeper reason for it. In our opinion the evidence 
given above that also the surface exponents extrapolate 
towards PC values in an intermediate regime strongly 
suggests that there is a genuine nonasymptotic PC-like 
regime, with a crossover to DP behavior at longer time 
scales. 

A prototype model in the PC class is the branching 
annihilating random walk with even offsprings (B ARWe) , 
defined by the reactions A 3A, 2A — > plus diffusion 
P, 0], which differs from the PCPD only for the reac- 
tion which creates particles. We argue that the early 
stages of the critical dynamics, when the system has a 
rather high particle density, are dominated by the an- 
nihilation process 2 A ^ 0, so that the substitution of 
the BARWe reaction A ^ 5A with that of the PCPD 
2A —f 3A may result in a very weak perturbation of the 
system. Therefore a transient PC-like regime may be 
observed for < ^ Tc, where Tc is some crossover time. 
This argument may help to explain features observed in 
the PCPD, and should be equally valid for other models 
where the annihilation is of the type 2A and with 
different creation rules nA {n+k)A with n > 2, fc > 0; 
for such systems we expect a transient PC regime as well. 

The study of reaction-diffusion systems where the an- 
nihilation and creation reactions involve n > 2 particles 
has recently drawn some attention |^ 0, |2^ . In partic- 
ular, we mention here the two cases recently considered 
by Odor [13 (i) 3A 5A, 2A ^ and (ii) AA 5A, 
2A 0. In model (i) he estimates ~ 0.28 (consis- 
tent with PC) for small diffusion rates and /3/i^|| ~ 0.24 
at stronger diffusion. Invoking some logarithmic correc- 
tions he claims that all values extrapolate to P/i'w ~ 0.22 
. In case (ii) the estimate is P/i^w ~ 0.28 both at high 
and low diffusions again consistent with the PC class 
value. The above observations suggest that these types 
of systems follow closely a critical behavior as described 
here for the PCPD, and it is thus plausible that they 
fall for sufficiently long times into the DP class. How- 
ever, it may turn out to be quite difficult to show this 
numerically, as we expect that increasing the number of 
particles involved in the creation and annihilation reac- 
tions will lead to models even harder to simulate and 
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FIG. 9: Plot of the surface effective exponent vs l/L for 
the parity conserving process A — > 3^4 and 2A — > in the case 
of reflecting boundary conditions. Inset: 7|ff vsl/L for ab- 
sorbing boundary conditions. Dashed lines are extrapolated 
curves through the DMRG data. Extrapolated values are in 
good agreement with Monte Carlo simulations results of Ref. 
[22 | (see text). 



reactions 2A A ^ 'HA and with single particle 

diffusion. We show how DMRG produces accurate sur- 
face critical exponents for this model, which are in good 
agreement with former Monte Carlo simulation results. 
For the single particle diffusion and pair annihilation we 
used the same rates as in Eq. while we assign a 

rate (1 - p){l - d) to the reaction 0^0 AAA. We re- 
strict ourselves to a single value of the diffusion constant 
d = 0.5. 

We first estimated the critical point aX p — pc ~ 
0.577(2) by means of Monte Carlo simulations using a 
standard approach. As mentioned above, for surface uni- 
versality in PC processes there are two possible types of 
boundary conditions leading to two distinct surface expo- 
nents j22| . In the first case, the system is truncated at one 
edge and no particles are allowed to cross the boundary 
site; we refer to this as reflecting boundary conditions. 
In the second case, particles are allowed to drop from the 
boundary. We implemented this type of boundary con- 
dition adding the boundary reaction A — s- (with rate d) 
which mimics the diffusion of particles out of the system. 
We refer to this implementation as absorbing boundary 
conditions. 



analyze than the PCPD. 

Very recently Kockelkoren and Chate analyzed a simi- 
lar set of models j^. In their formulation the fermionic 
constraint of only one particle per site is released. All 
the reactions of the type nA {n + k)A and 2 A 
with n > 2 were found to belong to the DP class. Sur- 
prisingly, in all those models the convergence to DP ex- 
ponents seems to be quite fast (at least for /3/v\^) and not 
plagued by the strong corrections found in the fermionic 
models. It would be interesting to study the same mod- 
els at different values of the diffusion constant, as in the 
PCPD the onset of crossover behavior is quite strongly 
influenced by the value of d. 

We are grateful to J.D. Noh, H. Hinrichsen, M. den 
Nijs, and F. van Wijland for useful discussions. 

APPENDIX: SURFACE CRITICAL BEHAVIOR 
IN THE PARITY CONSERVING PROCESS 2A ^ 
0, A ^ 3A 

We present here some results on the surface critical 
behavior of the parity conserving process defined by the 



Figure 121 shows the effective surface exponent j^g ver- 
sus 1/L in the case of reflecting boundary conditions, cal- 
culated both from the particle (circles) and pair (squares) 
densities. The same quantities are plotted in the inset in 
the case of absorbing boundary conditions. Notice that 
indeed the results confirm the existence of two distinct 
sets of surface exponents and that the data from pairs 
and particles merge for sufficiently long chains, indicat- 
ing that both quantities decay with the same exponent. 
Our estimates (3'^ /v± ~ 0.720(2) in the former case and 
(3^ /h'i_ ~ 1.10(1) in the latter are obtained from a poly- 
nomial extrapolation in 1/L. As the finite-size effects are 
rather small (see Fig. O , the extrapolated values are not 
very sensitive to the type of correction to scaling term 
used in the extrapolation. 

The Monte Carlo simulation results for the critical 
exponents are = 1.34(2) and = 2.04(2), for inactive 
and active boundary conditions, respectively. Combining 
these results with the PC class correlation length expo- 
nent = 1.83(3) Q, one finds P^/v^ = 0.73(1) (reflect- 
ing BCs) /3^/i^± = 1.11(1) (absorbing BCs), in very good 
agreement with the DMRG calculations. 
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